Tarea de Modelo lineales y de Sobrevivencia
Tarea 1

Parte 1

Librerias:

library(dplyr)
library(psych)
library(GGally)
library(knitr)
library(ggplot2)
library(gridExtra)
library(lmtest)
library(corrplot)
library(car)
library(datasets)
library(plotly)
library(nortest)
library(tseries)
library(carData)
library(MASS)

Del siguiente link: https://rpubs.com/Joaquin_AR/226291 ,realice las siguientes acciones:

a. Replique en un Script de R, los Ejemplos 1 y 2, usando todas las funciones que se presenta el artículo.

Ejemplo 1:

Para este primer ejemplo se desea generar un modelo que permita predecir la esperanza de vida media de los habitantes de una ciudad en función de diferentes variables. Se dispone de información sobre: habitantes, analfabetismo, ingresos, esperanza de vida, asesinatos, universitarios, heladas, área y densidad poblacional. Los datos se cargan acontinuación:

datos <- as.data.frame(state.x77)
datos <- rename(habitantes = Population, analfabetismo = Illiteracy,
                ingresos = Income, esp_vida = `Life Exp`, asesinatos = Murder,
                universitarios = `HS Grad`, heladas = Frost, area = Area,
                .data = datos)
datos <- mutate(.data = datos, densidad_pobl = habitantes * 1000 / area)

1) Analizar la relación entre variables.

Acontinucación, se estudia la relación entre las variables de la base de datos, esto para determinar si existe colinealidad o si se presentan varibales con relaciones de tipo no lineal. Para esto se corre la función cor() que devuelve las correlaciones entre variables y se redondea el resultado en el tercer decimal.

as.data.frame(round(cor(x = datos, method = "pearson"), 3))
##                habitantes ingresos analfabetismo esp_vida asesinatos
## habitantes          1.000    0.208         0.108   -0.068      0.344
## ingresos            0.208    1.000        -0.437    0.340     -0.230
## analfabetismo       0.108   -0.437         1.000   -0.588      0.703
## esp_vida           -0.068    0.340        -0.588    1.000     -0.781
## asesinatos          0.344   -0.230         0.703   -0.781      1.000
## universitarios     -0.098    0.620        -0.657    0.582     -0.488
## heladas            -0.332    0.226        -0.672    0.262     -0.539
## area                0.023    0.363         0.077   -0.107      0.228
## densidad_pobl       0.246    0.330         0.009    0.091     -0.185
##                universitarios heladas   area densidad_pobl
## habitantes             -0.098  -0.332  0.023         0.246
## ingresos                0.620   0.226  0.363         0.330
## analfabetismo          -0.657  -0.672  0.077         0.009
## esp_vida                0.582   0.262 -0.107         0.091
## asesinatos             -0.488  -0.539  0.228        -0.185
## universitarios          1.000   0.367  0.334        -0.088
## heladas                 0.367   1.000  0.059         0.002
## area                    0.334   0.059  1.000        -0.341
## densidad_pobl          -0.088   0.002 -0.341         1.000

Además, se representa la distribución de las vaiables mediante histogramas. Esto con ayuda de la fucnción multi.hist().

multi.hist(x = datos, dcol = c("blue", "red"), dlty = c("dotted", "solid"), global = F)

Otro ejemplo en el que se calculan las correlaciones y se generan histogramas a la vez, se encuentra con la función ggpairs().

ggpairs(datos, lower = list(continuous = "smooth"),
        diag = list(continuous = "barDiag"), axisLabels = "none")

De lo anterior se extra que las variables con mayor relación con esp_vida son asesinatos, analfabetismo y universitarios. Además, dado que las primeras dos variables poseen cierto grado de correlación, pude suceder que no sea útil intoducir ambas al modelo.

2) Generar el modelo.

Ahora, que se conoce la información repecto a las correlaciones y posibles distribuciones, se procede a generar un modelo mediante el método mixto y utilizando la medición AIC para la validación de los predictores.

modelo <- lm(esp_vida ~ habitantes + ingresos + analfabetismo + asesinatos +
             universitarios + heladas + area + densidad_pobl, data = datos )

summary(modelo)
## 
## Call:
## lm(formula = esp_vida ~ habitantes + ingresos + analfabetismo + 
##     asesinatos + universitarios + heladas + area + densidad_pobl, 
##     data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.47514 -0.45887 -0.06352  0.59362  1.21823 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     6.995e+01  1.843e+00  37.956  < 2e-16 ***
## habitantes      6.480e-05  3.001e-05   2.159   0.0367 *  
## ingresos        2.701e-04  3.087e-04   0.875   0.3867    
## analfabetismo   3.029e-01  4.024e-01   0.753   0.4559    
## asesinatos     -3.286e-01  4.941e-02  -6.652 5.12e-08 ***
## universitarios  4.291e-02  2.332e-02   1.840   0.0730 .  
## heladas        -4.580e-03  3.189e-03  -1.436   0.1585    
## area           -1.558e-06  1.914e-06  -0.814   0.4205    
## densidad_pobl  -1.105e-03  7.312e-04  -1.511   0.1385    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7337 on 41 degrees of freedom
## Multiple R-squared:  0.7501, Adjusted R-squared:  0.7013 
## F-statistic: 15.38 on 8 and 41 DF,  p-value: 3.787e-10

Con lo anterior se obtiene que el modelo pude explicar aproximadamente el 75% de la variabilidad, esto ya que el \(R^{2}\) fue de 0.7501. Por otro lado, ya que el p-valor obtenido fue significativo se pude aceptar que no se llegó al modelo por azar. Ademán, al menos uno de los coeficientes parciales de regresión es diferente de cero.

3) Selección de los mejores predictores.

Ahora se procede a seleccionar los mejores predictores empleando el método AIC.

step(object = modelo, direction = "both", trace = 1)
## Start:  AIC=-22.89
## esp_vida ~ habitantes + ingresos + analfabetismo + asesinatos + 
##     universitarios + heladas + area + densidad_pobl
## 
##                  Df Sum of Sq    RSS     AIC
## - analfabetismo   1    0.3050 22.373 -24.208
## - area            1    0.3564 22.425 -24.093
## - ingresos        1    0.4120 22.480 -23.969
## <none>                        22.068 -22.894
## - heladas         1    1.1102 23.178 -22.440
## - densidad_pobl   1    1.2288 23.297 -22.185
## - universitarios  1    1.8225 23.891 -20.926
## - habitantes      1    2.5095 24.578 -19.509
## - asesinatos      1   23.8173 45.886  11.707
## 
## Step:  AIC=-24.21
## esp_vida ~ habitantes + ingresos + asesinatos + universitarios + 
##     heladas + area + densidad_pobl
## 
##                  Df Sum of Sq    RSS     AIC
## - area            1    0.1427 22.516 -25.890
## - ingresos        1    0.2316 22.605 -25.693
## <none>                        22.373 -24.208
## - densidad_pobl   1    0.9286 23.302 -24.174
## - universitarios  1    1.5218 23.895 -22.918
## + analfabetismo   1    0.3050 22.068 -22.894
## - habitantes      1    2.2047 24.578 -21.509
## - heladas         1    3.1324 25.506 -19.656
## - asesinatos      1   26.7071 49.080  13.072
## 
## Step:  AIC=-25.89
## esp_vida ~ habitantes + ingresos + asesinatos + universitarios + 
##     heladas + densidad_pobl
## 
##                  Df Sum of Sq    RSS     AIC
## - ingresos        1     0.132 22.648 -27.598
## - densidad_pobl   1     0.786 23.302 -26.174
## <none>                        22.516 -25.890
## - universitarios  1     1.424 23.940 -24.824
## + area            1     0.143 22.373 -24.208
## + analfabetismo   1     0.091 22.425 -24.093
## - habitantes      1     2.332 24.848 -22.962
## - heladas         1     3.304 25.820 -21.043
## - asesinatos      1    32.779 55.295  17.033
## 
## Step:  AIC=-27.6
## esp_vida ~ habitantes + asesinatos + universitarios + heladas + 
##     densidad_pobl
## 
##                  Df Sum of Sq    RSS     AIC
## - densidad_pobl   1     0.660 23.308 -28.161
## <none>                        22.648 -27.598
## + ingresos        1     0.132 22.516 -25.890
## + analfabetismo   1     0.061 22.587 -25.732
## + area            1     0.043 22.605 -25.693
## - habitantes      1     2.659 25.307 -24.046
## - heladas         1     3.179 25.827 -23.030
## - universitarios  1     3.966 26.614 -21.529
## - asesinatos      1    33.626 56.274  15.910
## 
## Step:  AIC=-28.16
## esp_vida ~ habitantes + asesinatos + universitarios + heladas
## 
##                  Df Sum of Sq    RSS     AIC
## <none>                        23.308 -28.161
## + densidad_pobl   1     0.660 22.648 -27.598
## + ingresos        1     0.006 23.302 -26.174
## + analfabetismo   1     0.004 23.304 -26.170
## + area            1     0.001 23.307 -26.163
## - habitantes      1     2.064 25.372 -25.920
## - heladas         1     3.122 26.430 -23.877
## - universitarios  1     5.112 28.420 -20.246
## - asesinatos      1    34.816 58.124  15.528
## 
## Call:
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios + 
##     heladas, data = datos)
## 
## Coefficients:
##    (Intercept)      habitantes      asesinatos  universitarios         heladas  
##      7.103e+01       5.014e-05      -3.001e-01       4.658e-02      -5.943e-03

Por lo tanto el mejor modelo resulta de tomar habitantes, asesinatos, universitarios y heladas como predictores.

modelo <- (lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
              heladas, data = datos))
summary(modelo)
## 
## Call:
## lm(formula = esp_vida ~ habitantes + asesinatos + universitarios + 
##     heladas, data = datos)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.47095 -0.53464 -0.03701  0.57621  1.50683 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     7.103e+01  9.529e-01  74.542  < 2e-16 ***
## habitantes      5.014e-05  2.512e-05   1.996  0.05201 .  
## asesinatos     -3.001e-01  3.661e-02  -8.199 1.77e-10 ***
## universitarios  4.658e-02  1.483e-02   3.142  0.00297 ** 
## heladas        -5.943e-03  2.421e-03  -2.455  0.01802 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7197 on 45 degrees of freedom
## Multiple R-squared:  0.736,  Adjusted R-squared:  0.7126 
## F-statistic: 31.37 on 4 and 45 DF,  p-value: 1.696e-12

Se muestra un intervalo de confianza para cada coeficiente parcial de regresión.

confint(lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +
            heladas, data = datos))
##                        2.5 %        97.5 %
## (Intercept)     6.910798e+01 72.9462729104
## habitantes     -4.543308e-07  0.0001007343
## asesinatos     -3.738840e-01 -0.2264135705
## universitarios  1.671901e-02  0.0764454870
## heladas        -1.081918e-02 -0.0010673977

4) Validación de condiciones para la regresión múltiple lineal.

Ahora se valida la condición de variabilidad constante para los residuos es decir la homocesdasticidad.

plot1 <- ggplot(data = datos, aes(habitantes, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()
plot2 <- ggplot(data = datos, aes(asesinatos, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()
plot3 <- ggplot(data = datos, aes(universitarios, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()
plot4 <- ggplot(data = datos, aes(heladas, modelo$residuals)) +
    geom_point() + geom_smooth(color = "firebrick") + geom_hline(yintercept = 0) +
    theme_bw()
grid.arrange(plot1, plot2, plot3, plot4)

Ahora se testea la distribución normal de los residuos.

qqnorm(modelo$residuals)
qqline(modelo$residuals)

shapiro.test(modelo$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo$residuals
## W = 0.97935, p-value = 0.525

Se concluye que el análisis gráfico y la prueba, confirman la normalidad de los residuos.

Cuando graficamos los residuos en función de los valores ajustados por el modelo, esperamos ver que los residuos se distribuyan aleatoriamente alrededor de cero, sin mostrar ningún patrón. Además, esperamos que la variabilidad de los residuos sea más o menos constante a lo largo del eje X. Sin embargo, si identificamos algún patrón específico en la distribución de los residuos, como una forma cónica o una mayor dispersión en los extremos, esto indica que la variabilidad de los residuos está relacionada con los valores ajustados por el modelo, lo que sugiere una falta de homocedasticidad en el mismo.

ggplot(data = datos, aes(modelo$fitted.values, modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick", se = FALSE) +
geom_hline(yintercept = 0) +
theme_bw()

bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 6.2721, df = 4, p-value = 0.1797

Por lo tanto no hay evidencia de ausencia de homocedasticidad.

Relación entre los predictores.

corrplot(cor(dplyr::select(datos, habitantes, asesinatos,universitarios,heladas)),
         method = "number", tl.col = "black")

Análisis de Inflación de la Varianza.

vif(modelo)
##     habitantes     asesinatos universitarios        heladas 
##       1.189835       1.727844       1.356791       1.498077

El resultado anterior señala que no hay evidencia de una correlación lineal muy alta ni inflación de varianza.

Autocorrelación.

dwt(modelo, alternative = "two.sided")
##  lag Autocorrelation D-W Statistic p-value
##    1      0.02867262      1.913997   0.818
##  Alternative hypothesis: rho != 0

El resultado anterior muestra que no existe evidencia de auto correlación.

5) Identificación de posibles valores atípicos o influyentes.

Identificación de posibles valores atípicos.

datos$studentized_residual <- rstudent(modelo)

ggplot(data = datos, aes(x = predict(modelo), y = abs(studentized_residual))) +
  geom_hline(yintercept = 3, color = "grey", linetype = "dashed") +
  
# se identifican en rojo observaciones con residuos estandarizados absolutos > 3
geom_point(aes(color = ifelse(abs(studentized_residual) > 3, 'red', 'black'))) +
  scale_color_identity() +
  labs(title = "Distribución de los residuos studentized",
       x = "Predicción modelo") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5))

which(abs(datos$studentized_residual) > 3)
## integer(0)

No se encuentran valores atípicos.

summary(influence.measures(modelo))
## Potentially influential observations of
##   lm(formula = esp_vida ~ habitantes + asesinatos + universitarios +      heladas, data = datos) :
## 
##            dfb.1_ dfb.hbtn dfb.assn dfb.unvr dfb.hlds dffit   cov.r   cook.d
## Alaska      0.41   0.18    -0.40    -0.35    -0.16    -0.50    1.36_*  0.05 
## California  0.04  -0.09     0.00    -0.04     0.03    -0.12    1.81_*  0.00 
## Hawaii     -0.03  -0.57    -0.28     0.66    -1.24_*   1.43_*  0.74    0.36 
## Nevada      0.40   0.14    -0.42    -0.29    -0.28    -0.52    1.46_*  0.05 
## New York    0.01  -0.06     0.00     0.00    -0.01    -0.07    1.44_*  0.00 
##            hat    
## Alaska      0.25  
## California  0.38_*
## Hawaii      0.24  
## Nevada      0.29  
## New York    0.23

En la tabla anterior se listan las observaciones que son significativamente influyentes en al menos uno de los predictores.

influencePlot(modelo)

##               StudRes        Hat       CookD
## California -0.1500614 0.38475924 0.002879053
## Hawaii      2.5430162 0.23979244 0.363778638
## Maine      -2.2012995 0.06424817 0.061301962
## Nevada     -0.8120831 0.28860921 0.053917754
## Washington -1.4895722 0.17168830 0.089555784

El análisis anterior muestra valores preocupantes para los valores Leverages o Distancia Cook (California y Maine). Un estudio más exhaustivo recomendaría rehacer el modelo sin dichas observaciones y comparar los resultados.

6) Conclusión.

En conclusión el modelo de regresión lineal múltiple sería: Esperanza de vida \(=5.014e^{-05}\)habitantes \(- 3.001e^{−01}\)asesinatos \(+4.658e^{−02}\)universitarios \(−5.943e^{−03}\)heladas. que es capaz de explicar el 73.6% de la variabilidad observada en la esperanza de vida.

Ejemplo 2

Este ejmplo se ejecuta mediante una base de datos con información de 30 libros. Se conoce del peso total de cada libro, el volumen que tiene y el tipo de tapas (duras o blandas). Se quiere generar un modelo lineal múltiple que permita predecir el peso de un libro en función de su volumen y del tipo de tapas.

datos <- data.frame(peso = c(800, 950, 1050, 350, 750, 600, 1075, 250, 700,
                             650, 975, 350, 950, 425, 725),
                    volumen = c(885, 1016, 1125, 239, 701, 641, 1228, 412, 953,
                                929, 1492, 419, 1010, 595, 1034),
                    tipo_tapas = c("duras", "duras", "duras", "duras", "duras", 
                                   "duras", "duras", "blandas", "blandas",
                                   "blandas", "blandas", "blandas", "blandas",
                                   "blandas", "blandas"))
head(datos, 4)
##   peso volumen tipo_tapas
## 1  800     885      duras
## 2  950    1016      duras
## 3 1050    1125      duras
## 4  350     239      duras

1) Primero se analiza las correlaciones entre la variable cualitativas.

datos$tipo_tapas <- as.factor(datos$tipo_tapas)
pairs(x = datos)

cor.test(datos$peso, datos$volumen, method = "pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  datos$peso and datos$volumen
## t = 7.271, df = 13, p-value = 6.262e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7090393 0.9651979
## sample estimates:
##       cor 
## 0.8958988
ggplot(data = datos, mapping=aes(x = tipo_tapas, y = peso, color=tipo_tapas)) +
geom_boxplot() +
geom_jitter(width = 0.1) +
theme_bw() + theme(legend.position = "none")

Con los análisis anteriores se muestra que existe una relación lineal entre la variable peso y volumen. Además la variable tipo_tapas aparenta influir en el peso.

2) Modelo lineal multiple.

modelo <- lm(peso ~ volumen + tipo_tapas, data = datos)
summary(modelo)
## 
## Call:
## lm(formula = peso ~ volumen + tipo_tapas, data = datos)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -110.10  -32.32  -16.10   28.93  210.95 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      13.91557   59.45408   0.234 0.818887    
## volumen           0.71795    0.06153  11.669  6.6e-08 ***
## tipo_tapasduras 184.04727   40.49420   4.545 0.000672 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 78.2 on 12 degrees of freedom
## Multiple R-squared:  0.9275, Adjusted R-squared:  0.9154 
## F-statistic: 76.73 on 2 and 12 DF,  p-value: 1.455e-07
confint(modelo)
##                        2.5 %      97.5 %
## (Intercept)     -115.6237330 143.4548774
## volumen            0.5839023   0.8520052
## tipo_tapasduras   95.8179902 272.2765525

Por lo anterior, el modelo es capaz de explicar aproximadamente el 92.7% de la variabilidad observada, en el peso. Además, por el F teste, se puede decir que la muestra fue significativa.

3) Elección de los predictores.

En este caso solo hay dos posibles predicotres que mediante la función summary() se constataron como importantes.

4) Condiciones para la regresión lineal múltiple.

Relación lineal entre los predictores numéricos y la variable dependiente:

ggplot(data = datos, aes(x = volumen, y = modelo$residuals)) +
geom_point() +
geom_smooth(color = "firebrick") +
geom_hline(yintercept = 0) +
theme_bw()

Se satisface la condición de linealidad y se observa un dato que podría ser atípico.

Distribución normal de los residuos:

qqnorm(modelo$residuals)
qqline(modelo$residuals)

shapiro.test(modelo$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo$residuals
## W = 0.85497, p-value = 0.02043

Por lo anterior, se concluye que la condición de normalidad no se cumple, posiblemente por el dato atípico observado anteriormente. Por lo tanto, se vuelve a correr el código sin el valor atípico.

dato_at <- which.max(modelo$residuals)
shapiro.test(modelo$residuals[-dato_at])
## 
##  Shapiro-Wilk normality test
## 
## data:  modelo$residuals[-dato_at]
## W = 0.9602, p-value = 0.7263

Por lo tanto, los residuos sí se distribuyen de forma normal.

Variabilidad constante de los residuos:

ggplot(data = data.frame(predict_values = predict(modelo),
                         residuos = residuals(modelo)),
       aes(x = predict_values, y = residuos)) +
    geom_point() +
    geom_smooth(color = "firebrick", se = FALSE) +
    geom_hline(yintercept = 0) +
    theme_bw()

bptest(modelo)
## 
##  studentized Breusch-Pagan test
## 
## data:  modelo
## BP = 2.0962, df = 2, p-value = 0.3506

No hay evidencias de falta de homocedasticidad. Además, dado que solo hay un predictor cuantitativo, no se puede dar colinealidad.

Autocorrelación:

dwt(modelo,alternative = "two.sided")
##  lag Autocorrelation D-W Statistic p-value
##    1    0.0004221711      1.970663   0.744
##  Alternative hypothesis: rho != 0

No hay evidencias de autocorrelación.

5) Identificación de posibles valores atípicos o influyentes

outlierTest(modelo)
##    rstudent unadjusted p-value Bonferroni p
## 13 5.126833         0.00032993     0.004949

Justo como se vio anteriormente, existe un dato atípico.

summary(influence.measures(modelo))
## Potentially influential observations of
##   lm(formula = peso ~ volumen + tipo_tapas, data = datos) :
## 
##    dfb.1_ dfb.vlmn dfb.tp_t dffit   cov.r   cook.d hat  
## 4  -0.16   0.18    -0.10    -0.23    1.98_*  0.02   0.36
## 11  0.70  -1.26_*   0.57    -1.54_*  0.83    0.64   0.38
## 13  0.31   0.67    -1.31_*   2.07_*  0.04_*  0.46   0.14
influencePlot(modelo)

##       StudRes       Hat      CookD
## 4  -0.3008849 0.3616894 0.01850173
## 11 -1.9897106 0.3757842 0.63729788
## 13  5.1268330 0.1397761 0.45819722

El análisis anterior muestra la existencia de múltiples observaciones influyentes, pero ninguna preocupante bajo los valores leverage hat o distancia Cook.

6) Conclusión.

El modelo Peso \(=13.91557+0.71795\)volumen \(+184.04727\)tipotapas, es capaz de explicar el 92.7% de la variabilidad observada en el pseo de los libros.

b.Explique brevemente las siguientes funciones, incluyendo valores recibe y cuáles son sus salidas.

  • multi.hist de la biblioteca psych: Recibe una matriz o dataframe y devuelve un plot con multiples histogramas para cada variable. Además, puede representar distribuciones de densidad para cada gráfico.
  • ggpairs de la biblioteca GGally: Esta función recibe un dataframe con variables continuas o categoricas y devuleve un gráfico que contiene las correlaciones encima de la diagonal, por debajo los diagramas de dispersión entre las variables continuas, en la diagonal los gráficos de densidad de las variables continuas y en los lados los histogramas y box plots de la combinación entre variables continuas y categóricas.
  • Función Step de la sección 3 del ejemplo 1: Esta función permite encontrar el mejor modelo basado en AIC utilizando a elección el método forward, backward o moxto. Esta función recibe un modelo el cual es un objeto de tipo lm o glm, y devuelve el modelo seleccionado.
  • Confint : Recibe un objeto de tipo modelo y devuelve una matriz o vector con los limites superior e inferior de confianza para cada parámetro.
  • Corrplot de la biblioteca corrplot: Esta función recibe una matriz de correlaciones y devuelve un gráfico de dichas correlaciones.
  • Dwt de la biblioteca car: Esta función recibe una serie de tiempo univariada o multivariada (vectores numéricos, matrices y dataframes son también aceptados) y calcula los coeficientes de la transformada wavelet discretaDevuelve. Devuelve un objeto de tipo dwt con información respecto a la autocorrelación.
  • InfluencePlot: Esta función recibe un modelo lineal y retorna un gráfico de “burbujas” de los residuos de Studentized frente a los valores hat, con las áreas de los círculos que representan las observaciones proporcionales al valor Distancia de Cook.
  • outlierTest de la biblioteca car: Esta función recibe un modelo de regresión lineal e informa de los valores atípicos basado en los p-valores de Bonferroni.
  • influence.measure: Esta función se utiliza para calcular algunos de los diagnósticos de regresión para modelos lineales y lineales generalizados. Recibe un modelo lineal y retorna objeto que al pasarlo por la función summary, revela información relevante sobre la influencia de las observaciones.

c.Explique los test incluyendo cual es la hipótesis nula.

  • Shapiro-Wilk para normalidad: Esta prueba estadística se utiliza para determinar si un conjunto de datos proviene de una distribución normal, para ello se establecen las hipótesis: \(H_{0}:\) La distribución es normal y \(H_{1}:\) La distribución no es normal. Al realizar la prueba y dado un nivel de significancia, el p-valor arrojado indicará si se rechaza la hipotesis nula (\(H_{0}\)) o sí, por lo contrario, se acepta.

  • Studentized Breusch-Pagan para homocedasticidad: Esta prueba se implementa para corroborar la homocedasticidad al analizar si la varianza estimada de los residuos de una regresión dependen de los valores de las variables independientes. En este test las hipótesis son: \(H_{0}:\) Los errores tiene varianza constante y \(H_{1}:\) Los errores no tiene varianza constante. Además, como en el test anterior dado un nivel de significancia, el p-valor arrojado podrá o no rechazar la hipótesis nula.

Parte 2

state.x77=as.data.frame(state.x77)
head(state.x77)
##            Population Income Illiteracy Life Exp Murder HS Grad Frost   Area
## Alabama          3615   3624        2.1    69.05   15.1    41.3    20  50708
## Alaska            365   6315        1.5    69.31   11.3    66.7   152 566432
## Arizona          2212   4530        1.8    70.55    7.8    58.1    15 113417
## Arkansas         2110   3378        1.9    70.66   10.1    39.9    65  51945
## California      21198   5114        1.1    71.71   10.3    62.6    20 156361
## Colorado         2541   4884        0.7    72.06    6.8    63.9   166 103766

1.Cálcule la matriz de Correlaciones de la base de datos states.x77

matriz_cor <- cor(state.x77)
matriz_cor
##             Population     Income  Illiteracy    Life Exp     Murder
## Population  1.00000000  0.2082276  0.10762237 -0.06805195  0.3436428
## Income      0.20822756  1.0000000 -0.43707519  0.34025534 -0.2300776
## Illiteracy  0.10762237 -0.4370752  1.00000000 -0.58847793  0.7029752
## Life Exp   -0.06805195  0.3402553 -0.58847793  1.00000000 -0.7808458
## Murder      0.34364275 -0.2300776  0.70297520 -0.78084575  1.0000000
## HS Grad    -0.09848975  0.6199323 -0.65718861  0.58221620 -0.4879710
## Frost      -0.33215245  0.2262822 -0.67194697  0.26206801 -0.5388834
## Area        0.02254384  0.3633154  0.07726113 -0.10733194  0.2283902
##                HS Grad      Frost        Area
## Population -0.09848975 -0.3321525  0.02254384
## Income      0.61993232  0.2262822  0.36331544
## Illiteracy -0.65718861 -0.6719470  0.07726113
## Life Exp    0.58221620  0.2620680 -0.10733194
## Murder     -0.48797102 -0.5388834  0.22839021
## HS Grad     1.00000000  0.3667797  0.33354187
## Frost       0.36677970  1.0000000  0.05922910
## Area        0.33354187  0.0592291  1.00000000

2. Genere tres modelos lineales usando la función lm de R para:

Modelo 1: Murder ~ Population

model1 <- lm(Murder ~ Population, data = state.x77)
summary(model1)
## 
## Call:
## lm(formula = Murder ~ Population, data = state.x77)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.9855 -3.0119 -0.3128  2.4986  7.9014 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.1713934  0.6869410   8.984 7.49e-12 ***
## Population  0.0002841  0.0001121   2.535   0.0146 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.503 on 48 degrees of freedom
## Multiple R-squared:  0.1181, Adjusted R-squared:  0.09972 
## F-statistic: 6.427 on 1 and 48 DF,  p-value: 0.01455

Modelo 2: Murder ~ Illiteracy

model2 <- lm(Murder ~ Illiteracy, data = state.x77)
summary(model2)
## 
## Call:
## lm(formula = Murder ~ Illiteracy, data = state.x77)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5315 -2.0602 -0.2503  1.6916  6.9745 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.3968     0.8184   2.928   0.0052 ** 
## Illiteracy    4.2575     0.6217   6.848 1.26e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.653 on 48 degrees of freedom
## Multiple R-squared:  0.4942, Adjusted R-squared:  0.4836 
## F-statistic: 46.89 on 1 and 48 DF,  p-value: 1.258e-08

Modelo 3: Murder ~ Population + Illiteracy

model3 <- lm(Murder ~ Population + Illiteracy, data = state.x77)
summary(model3)
## 
## Call:
## lm(formula = Murder ~ Population + Illiteracy, data = state.x77)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7652 -1.6561 -0.0898  1.4570  7.6758 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.652e+00  8.101e-01   2.039  0.04713 *  
## Population  2.242e-04  7.984e-05   2.808  0.00724 ** 
## Illiteracy  4.081e+00  5.848e-01   6.978 8.83e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.481 on 47 degrees of freedom
## Multiple R-squared:  0.5668, Adjusted R-squared:  0.5484 
## F-statistic: 30.75 on 2 and 47 DF,  p-value: 2.893e-09

3. Represente gráficamente los gráficos de dispersión correspondientes, las rectas de mínimos cuadrados para los tres modelos anteriores

de mínimos cuadrados para los tres modelos anteriores.

# Gráfico de dispersión y recta de mínimos cuadrados para el modelo 1: Murder ~ Population
plot(state.x77$Population, state.x77$Murder, main = "Murder vs Population", xlab = "Population", ylab = "Murder")
abline(model1, col = "red")

# Gráfico de dispersión y recta de mínimos cuadrados para el modelo 2: Murder ~ Illiteracy
plot(state.x77$Illiteracy, state.x77$Murder, main = "Murder vs Illiteracy", xlab = "Illiteracy", ylab = "Murder")
abline(model2, col = "blue")

# Gráfico de dispersión y recta de mínimos cuadrados para el modelo 3: Murder ~ Population + Illiteracy

# Crear una cuadrícula de valores para Illiteracy y Population
illiteracy_grid <- seq(min(state.x77$Illiteracy), max(state.x77$Illiteracy), length.out = 20)
population_grid <- seq(min(state.x77$Population), max(state.x77$Population), length.out = 20)
grid <- expand.grid(Illiteracy = illiteracy_grid, Population = population_grid)

# Calcular los valores predichos de Murder para cada combinación de Illiteracy y Population
predictions <- predict(model3, newdata = grid)

# Calcular los valores predichos de Murder para los datos originales
predictions_original <- predict(model3)

# Coeficientes del plano de mínimos cuadrados
intercept <- coef(model3)[1]
slope_population <- coef(model3)[2]
slope_illiteracy <- coef(model3)[3]

# Calcular los valores de Murder para la recta de mínimos cuadrados
murder_line <- intercept + slope_population * grid$Population + slope_illiteracy * grid$Illiteracy

# Crear el gráfico 3D con plotly
scatter_plot <- plot_ly(x = grid$Illiteracy, y = grid$Population, z = predictions, type = "surface")
scatter_plot <- add_trace(scatter_plot, x = state.x77$Illiteracy, y = state.x77$Population, z = state.x77$Murder, mode = "markers", type = "scatter3d", marker = list(size = 5))

# Agregar plano de mínimos cuadrados
scatter_plot <- add_trace(scatter_plot, x = grid$Illiteracy, y = grid$Population, z = murder_line, mode = "lines", type = "scatter3d", line = list(color = "blue", width = 3))

# Etiquetas y título
scatter_plot <- layout(scatter_plot, scene = list(xaxis = list(title = "Illiteracy"), yaxis = list(title = "Population"), zaxis = list(title = "Murder")), title = "Murder vs Population + Illiteracy")

# Mostrar el gráfico
scatter_plot

4. Evalue para el modelo Murder ~ Population las hipótesis de homocedasticidad, normalidad de errores y ausencia de puntos influyentes o atípicos.

residuos <- residuals(model1)
plot(residuos)

Homocedasticidad

# Prueba Breusch-Pagan para la homocedasticidad
bptest(model1)
## 
##  studentized Breusch-Pagan test
## 
## data:  model1
## BP = 1.4303, df = 1, p-value = 0.2317

Como el valor p (0.2317) es mayor que el nivel de significancia típico de 0.05, no hay suficiente evidencia para rechazar la hipótesis nula de homocedasticidad.Por ende, según esta prueba, el modelo Murder ~ Population cumple con el supuesto de homocedasticidad.

Normalidad de errores

# Prueba Shapiro-Wilk para la Normalidad
shapiro.test(residuos)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos
## W = 0.95662, p-value = 0.0642

Dado que el valor p (0.0642) es mayor que el nivel de significancia típico de 0.05, no hay suficiente evidencia para rechazar la hipótesis nula de normalidad. Por lo tanto, según esta prueba, los residuos del modelo parecen seguir una distribución normal.

# Prueba Kolmogorov-Smirnov para la Normalidad
lillie.test(residuos)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  residuos
## D = 0.080377, p-value = 0.5801

Dado que el valor p (0.5801) es mayor que el nivel de significancia comúnmente utilizado de 0.05, no hay suficiente evidencia para rechazar la hipótesis nula de normalidad. Por lo tanto, según esta prueba, los residuos del modelo parecen seguir una distribución normal.

# Prueba Jarque Bera para la Normalidad
jarque.bera.test(residuos)
## 
##  Jarque Bera Test
## 
## data:  residuos
## X-squared = 2.603, df = 2, p-value = 0.2721

Dado que el valor p (0.2721) es mayor que el nivel de significancia típicamente utilizado de 0.05, no hay suficiente evidencia para rechazar la hipótesis nula de que los datos provienen de una distribución normal. Por lo tanto, en este caso, los datos podrían considerarse aproximadamente normales según el test de Jarque-Bera.

Puntos influyentes o atípicos

qqnorm(residuos)

num_bins <- 30
hist(residuos, breaks = num_bins, col = "skyblue", main = "Histograma de Residuos", xlab = "Residuos")

mediaRes <- mean(residuos)
desRes <- sd(residuos)

valorRefAtipD <- mediaRes+2*desRes
valorRefAtipI <- mediaRes-2*desRes

atipDerecha <- residuos[residuos > valorRefAtipD]
atipIzquierda <- residuos[residuos < valorRefAtipI]

# Imprime los residuos atípicos
print("Residuos atípicos en el extremo derecho:")
## [1] "Residuos atípicos en el extremo derecho:"
print(atipDerecha)
##  Alabama 
## 7.901416
print("Residuos atípicos en el extremo izquierdo:")
## [1] "Residuos atípicos en el extremo izquierdo:"
print(atipIzquierda)
## named numeric(0)
# Prueba de los residuos estandarizados
residuos_estandarizados <- residuos / sd(residuos)
valores_atipicos_residuos <- which(abs(residuos_estandarizados) > 2)  # Umbral común: 2 desviaciones estándar
print(valores_atipicos_residuos)
## Alabama 
##       1
outlier_test <- outlierTest(model1)
print(outlier_test)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
##         rstudent unadjusted p-value Bonferroni p
## Alabama 2.388285           0.020997           NA
plot(model1)

Según los gráficos y pruebas anteriores, hay evidencia de al menos un valor atípico (Alabama) por la derecha.

influencePlot(model1)

##               StudRes        Hat       CookD
## Alabama     2.3882848 0.02040822 0.054112891
## California -0.6492747 0.31422549 0.097758047
## Georgia     1.8723014 0.02047985 0.034828716
## New York   -0.1300681 0.21582941 0.002376835

Además según lo anteior, se presentan ciertos valores influyentes.

5. Estima las tasas de homicidio para niveles de Illiteracy de 0.25, 1.2, 2.1, 3.0, 4.0. en el segundo modelo.

coeficientes <- coef(model2)

niveles_illiteracy <- c(0.25, 1.2, 2.1, 3.0, 4.0)

tasas_homicidioE <- coeficientes[1] + coeficientes[2] * niveles_illiteracy

resultado <- data.frame(Illiteracy = niveles_illiteracy, Tasa_Homicidio_Estimada = tasas_homicidioE)
print(resultado )
##   Illiteracy Tasa_Homicidio_Estimada
## 1       0.25                3.461140
## 2       1.20                7.505724
## 3       2.10               11.337435
## 4       3.00               15.169146
## 5       4.00               19.426603

Parte 3

Se carga la base de datos correspondiente

Y = c(47.83,59.51,50.38,53.24,47.26,50.52)
Ynames = c("obs1","obs2","obs3","obs4","obs5","obs6")
names = c("b0", "tasaFlujoGas1", "tasaFlujoGas2", "aperturaBoq", "tempGas")

X = matrix(c(1, 124.17,134.07,23.32,210.11, 
    1,149.46,142.21,41.82,229.67,
    1,131.65,146.62,21.14,231.10,
    1,139.49,136.16,45.79,206.03,
    1,113.03,125.41,41.51,222.67,
    1,134.57,165.84,32.42,219.59),nrow=6,byrow = TRUE,dimnames=list(Ynames,names))

X2 = data.frame(cbind(Y,X[,-1]))

a. Estime un modelo de regresión lineal, con la función lm() para la variable X2.

regresion_lineal = lm(Y ~ ., data = X2)
summary(regresion_lineal)
## 
## Call:
## lm(formula = Y ~ ., data = X2)
## 
## Residuals:
##     obs1     obs2     obs3     obs4     obs5     obs6 
##  0.53326  0.48658 -0.67440 -0.57468  0.03905  0.19019 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -10.23386   12.59988  -0.812   0.5657  
## tasaFlujoGas1   0.33540    0.05226   6.417   0.0984 .
## tasaFlujoGas2  -0.07329    0.04757  -1.541   0.3665  
## aperturaBoq     0.08865    0.05935   1.494   0.3756  
## tempGas         0.11253    0.05326   2.113   0.2814  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.159 on 1 degrees of freedom
## Multiple R-squared:  0.9867, Adjusted R-squared:  0.9334 
## F-statistic: 18.51 on 4 and 1 DF,  p-value: 0.1724

b. La función lm() estima la inversa de XtX con base en el proceso de descomposición QR, de Gram–Schmidt, realice un breve descripción de este proceso y verifique usando la función qr.solve () que los estimadores del resultado de la regresión lineal de punto a son los mismos.

El proceso de descomposición QR se desarrolla de la siguiente forma:

Sea \(X\) una matriz nxm . Si consideramos las m columnas de \(X\), denotadas por \(x_j\) podemos obtener una base ortogonal de vectores unitarios \(q_j\) tales que:

Si \(X\) tiene rango completo (es decir, sus columnas son linealmente independientes), podemos relacionar los vectores \(x_j\) y \(q_j\) mediante la ecuación

\(x_k = r_{1k}q_1 + r_{2k}q_2 + r_{3k}q_3 + ... + r_{kk}q_k.\), con k = 1, …,m

Como los vectores \(q_j\) son ortogonales, los coeficientes vienen dados por

\(rij = q_i*a_j\) con \(i \not = j\)

Esto permite escribir

\(X = QR\) , donde \(Q\) es una matriz nxm que satisface que \(QQ = I\). Y, \(R\) es una matriz mxm triangular superior.

Si \(X\) no tiene rango completo, entonces basta con seleccionar \(q_k\) de manera tal que sea ortogonal a \((q_1, . . . , q_{k-1})\). De esta manera también se puede escribir \(X = QR\).

Finalmente, la inversa de la matriz \(X^tX\) se puede calcular mediante la descomposición \(QR\) de la matriz \(X\).

Verificación usando qr.solve()

Empleando la función qr.solve se obtiene lo siguiente:

QR = qr(X)
qr.solve(QR, Y)
##            b0 tasaFlujoGas1 tasaFlujoGas2   aperturaBoq       tempGas 
##  -10.23386366    0.33540014   -0.07329371    0.08864733    0.11252772

Como se puede observar, son los mismos valores de beta que brinda la función lm en el punto a.

c. Usando una matriz pseudo-inversa (Moonre-Penrose) y según lo visto en clases, reconstruya los valores de la regresión del punto a.

  1. Estimadores

    Primeramente, se calcula la matriz ortogonal

X_matriz = as.matrix(X) 
XTX = t(X_matriz) %*% X_matriz 
XTXI = ginv(XTX) 
XTXI
##              [,1]          [,2]          [,3]          [,4]          [,5]
## [1,] 118.13156142 -0.0340693673 -0.1000545541 -0.1742167735 -0.4243787993
## [2,]  -0.03406937  0.0020325502 -0.0010235978 -0.0010630811 -0.0002401096
## [3,]  -0.10005455 -0.0010235978  0.0016840371  0.0009228886 -0.0001597012
## [4,]  -0.17421677 -0.0010630811  0.0009228886  0.0026209996  0.0004267736
## [5,]  -0.42437880 -0.0002401096 -0.0001597012  0.0004267736  0.0021107277
La solución de los betas es la siguiente:
Betas = XTXI %*% t(X_matriz) %*% Y
Betas
##              [,1]
## [1,] -10.23386366
## [2,]   0.33540014
## [3,]  -0.07329371
## [4,]   0.08864733
## [5,]   0.11252772
Los cuales son los mismo que los obtenidos mediante la función lm.

  1. Std error de los estimadores

    Como primer paso, se calcula la matriz de covarianzas

Var_Betas = vcov(regresion_lineal) 
Var_Betas
##                (Intercept) tasaFlujoGas1 tasaFlujoGas2  aperturaBoq
## (Intercept)   158.75700508 -0.0457858226 -0.1344633150 -0.234129921
## tasaFlujoGas1  -0.04578582  0.0027315443 -0.0013756131 -0.001428675
## tasaFlujoGas2  -0.13446331 -0.0013756131  0.0022631774  0.001240270
## aperturaBoq    -0.23412992 -0.0014286747  0.0012402700  0.003522361
## tempGas        -0.57032267 -0.0003226833 -0.0002146225  0.000573541
##                     tempGas
## (Intercept)   -0.5703226673
## tasaFlujoGas1 -0.0003226833
## tasaFlujoGas2 -0.0002146225
## aperturaBoq    0.0005735410
## tempGas        0.0028366070
Posteriormente, se obtiene las std error de los betas
errores = c()

for (i in 1: 5) {
  errores[i] = sqrt(Var_Betas[i,i]) 
}  

errores_Betas = matrix(errores, dimnames = list(names, "Std error" )) 
errores_Betas
##                 Std error
## b0            12.59988115
## tasaFlujoGas1  0.05226418
## tasaFlujoGas2  0.04757286
## aperturaBoq    0.05934948
## tempGas        0.05325981
Como se puede observar, son iguales a los errores calculados con la
función lm.

  1. T-value

    Se procede a calcular el los T-valores de cada variable:

t_valor_matriz = c() 

for (i in 1: 5) {  
  t_valor_matriz[i] = Betas[i,]/errores_Betas[i,] 
} 

t_valores = matrix(t_valor_matriz,dimnames = list(names, "T-value" ))
t_valores
##                  T-value
## b0            -0.8122191
## tasaFlujoGas1  6.4174002
## tasaFlujoGas2 -1.5406622
## aperturaBoq    1.4936495
## tempGas        2.1128076

  1. Pr(>|t|)

    Con lo anterior, se obtiene que las probabilidades son:

    probabilidades = c() 

for (i in 1: 5) {
probabilidades[i] = 2*(1-pt(abs(t_valor_matriz[i]), 1)) 
} 

P_T = matrix(probabilidades,dimnames = list(names, "Pr(>|t|)" )) 
P_T
    ##                 Pr(>|t|)
## b0            0.56573154
## tasaFlujoGas1 0.09841069
## tasaFlujoGas2 0.36651617
## aperturaBoq   0.37558169
## tempGas       0.28142638

    De tal manera, comparando con lo obtenido por lm, se tienen los mismos resultados.

  2. Residual estándar error.

    Primero, se obtiene la matriz de proyección P

P = X %*% XTXI %*% t(X)
P
##             obs1        obs2        obs3        obs4         obs5         obs6
## obs1  0.78840211 -0.19307656  0.26760326  0.22803272 -0.015495472 -0.075466065
## obs2 -0.19307656  0.82382358  0.24417974  0.20807284 -0.014139141 -0.068860462
## obs3  0.26760326  0.24417974  0.66156796 -0.28838804  0.019596787  0.095440297
## obs4  0.22803272  0.20807284 -0.28838804  0.75425595  0.016699007  0.081327525
## obs5 -0.01549547 -0.01413914  0.01959679  0.01669901  0.998865255 -0.005526437
## obs6 -0.07546607 -0.06886046  0.09544030  0.08132752 -0.005526437  0.973085142
Luego, se calcula la matriz IP
Id = diag(rep(1,6)) 
IP = Id-P
Tercero, se obtiene la SSE (suma de los cuadrados residuales)
SSE=t(Y)%*%IP%*%Y
Finalmente, se calcula el error residual estándar
error_residual = sqrt(SSE/(6-5)) 
error_residual
##          [,1]
## [1,] 1.159267
La cual es la misma que da la función lm.

  1. R cuadrado

    Primeramente, se obtiene la SST (total de las sumas la cuadrado)

SST = t(Y)%*%(Y) 
SST
##          [,1]
## [1,] 15987.57
Se centran las observaciones, por lo que se tiene el SSTC (SST
centrado)
y_mean = mean(Y) 
SSTC = SST-6*y_mean^2 
SSTC
##          [,1]
## [1,] 100.8377
Por consiguiente, el R cuadrado es:
R2= 1-SSE/SSTC 
R2
##           [,1]
## [1,] 0.9866726
que es el mismo que el dado por la función lm.
  1. R cuadrado Ajustado
R2adj = 1-((6-1)/(6-5))*(1-R2) 
R2adj
##           [,1]
## [1,] 0.9333632
 Es posible identificar que es igual al calculado por la función lm.

  1. F-estadistico y su p-value.

    Para obtener el F-estadístico, es necesario conocer sus grados de libertad. Se sabe que posee t y N-r = 1 grados de libertad. Entonces, queda pendiente conocer t.

    Por medio de las probabilidades Pr(>|t|), se puede identificar que la variable tasaFlujoGas1 posee la menor probabilidad y por ende se considera como una variable significativa junto a B0, con respecto a las otras. Dado que t es la cantidad de variables que no aportan información significativa, se tiene que t = 4.

    De tal manera, el estadístico F es:

F = (R2/4)/(1-R2) 
F
##          [,1]
## [1,] 18.50841
  Finalmente, el p-valor corresponde a :
p_valor = 1-pf(F,4,1) 
p_valor
##           [,1]
## [1,] 0.1723969
  Como se puede ver,ambos valores son los mismos que los dados por
  la función lm.

d. Realice una prueba de contraste para determinar si se acepta o no la hipótesis que Tasa Flujo Gas 2 + aperturaBoq = 0.

Se realiza una prueba de contraste con las siguientes hipótesis:

  • \(H_0 : B2 + B3 = 0\)
  • \(H_1 : B2 + B3 \not = 0\)

Primero, se calcular el error:

error_suma =  sqrt(0.0022631774 + 0.003522361 +2*0.001240270)
error_suma
## [1] 0.09091798

Así, el T-valor es:

t_valor_suma = (-0.07329 + 0.08865)/error_suma
t_valor_suma
## [1] 0.1689435

Finalmente, se calcula el p-valor para determinar si se rechaza o no la hipótesis nula. Para este, tenemos r = 5 y N = 6. De tal manera, el p-valor es:

p_valor_suma = 2*(1-pt(t_valor_suma, 6-5))
p_valor_suma
## [1] 0.8934533

Por tanto, con un nivel de confianza del 95%, se acepta la hipótesis nula.